en 

o 

(N 



< 



< 



> 

O 



en 
en 

o 



c^ 



An explicit solution for implicit time stepping 
in finite strain viscoelasticity. 



A.V. Shutov*, R. Landgraf, J. Thiemann 



Chemnitz University of Technology, Department of Solid Mechanics, Str. d. 
Oj. Nationen 62, D- 09 111 Chemnitz, Cermany 



Abstract 

We consider the numerical treatment of one of the most popular finite strain 
models of the viscoelastic Maxwell body. This model is based on the multiplicative 
decomposition of the deformation gradient, combined with Neo-Hookean hypere- 
lastic relations between stresses and elastic strains. The evolution equation is six 
dimensional. For the corresponding local initial value problem, a fully implicit in- 
tegration procedure is considered, and a simple explicit update formula is derived. 
Thus, no local iterative procedure is required, which makes the numerical scheme 



00 ' more robust and efficient. The resulting integration algorithm is unconditionally 



stable and first order accurate. The incompressibility constraint of the inelastic 
flow is exactly preserved. A rigorous proof of the symmetry of the consistent tan- 
gent operator is provided. Moreover, some properties of the numerical solution, like 

ff^ ■ invariance under the change of the reference configuration and positive energy dis- 

sipation within a time step, are discussed. Numerical tests show that, in terms of 
accuracy, the proposed integration algorithm is equivalent to the classical implicit 

^ ' scheme based on the exponential mapping. Finally, in order to check the stability of 

^ . the algorithm numerically, a representative initial boundary value problem involving 

finite viscoelastic deformations is considered. An FEM solution of the representative 
problem using MSC.MARC is presented. 
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1 Introduction 



Among idealized models of linear viscoelasticity, the so-called Maxwell fluid (MF) is com- 
monly encountered in material modelling [53J.6,2]. The one- dimensional rheological inter- 
pretation of this model is shown in Fig. la[^ A series of Maxwell elements connected in 
parallel [67] can be utilized to represent viscoelastic properties of polymers (Fig. lb). In 
that case, the stresses acting in the Maxwell elements can be associated with overstresses 
[16]. Next, a slightly modifled Maxwell element can be adopted to capture the nonlinear 
kinematic hardening in metals (Fig. Ic). In that case, the corresponding Maxwell stresses 
are interpreted as backstresses [39,3,66,7,57]. Moreover, within some phenomenological 
approaches to metal plasticity, the distortional hardening in metals can be captured us- 
ing the modifled MF [58,60]. Other groups of materials like shape memory alloys [18,20] 
and biological tissues [8] can be modelled using MF. Some apphcations of MF to flnite 
deformations of geological structures [61,49] and to fluid mechanics [1] are reported in the 
literature as well. 
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Fig. 1. a) A one-dimensional Maxwell body consists of an elastic spring (Hooke body) coupled 
in series with a viscous dashpot (Newton body), b) Generalized Maxwell body, also known as 
Wiechert model (or Zener model in a special case) used for better description of viscoelastic 
properties, c) modified Schwedoff model used to represent nonlinear kinematic hardening. 

In the flnite strain range, numerous constitutive models of the MF exist (see, among 
others, [37,27,33,48,46,47,23,4,51,1,16,17]). Different variants were compared through nu- 
merical tests in [10,34]. In this paper we consider one of the most popular models of 
the MF. The corresponding constitutive equations are summarized in Section 2 of this 
work. The model under consideration is a special case of the flnite strain viscoplasticity 
model proposed by Simo and Miehe [62], and it has the same structure as the well known 
model of associative elastoplasticity considered by Simo [63]. These models were devel- 
oped within the framework of multiplicative inelasticity in combination with hyperelastic 
constitutive relations. The corresponding inelastic flow rule ^ I is six dimensional since the 
inelastic spin plays no role due to elastic isotropy. A version of the MF which is equivalent 
to the version of Simo & Miehe was considered later in material (Lagrangian) description 
by Lion [38]. This Lagrangian formulation was adopted in [40,18,12,54,56,58]. The spatial 



^ A two-dimensional rheological model of the Maxwell fluid and its modifications can be found 
in [58,60]. 

^ In this paper, the evolution equation is referred to as "inelastic flow rule" in order to stress 
that the model is a special case of a viscoplastic model. 



(Eulerian) constitutive equations proposed by Simo & Miehe were utilized later in the 
comprehensive study by Reese and Govindjee [51], and by many others (see, for instance, 
[24,43,49,30,15,22,50,36]). 

In the modern literature on numerics, much attention is paid to general procedures which 
can be implemented to different types of constitutive relations in a straightforward man- 
ner. Because of their generality, such procedures are not always efficient being compared 
to algorithms which make use of the special structure of the underlying constitutive equa- 
tions. Due to the high prevalence of the Simo & Miehe version of the Maxwell fluid in 
material modelling, efficient and robust numerical integration of the underlying evolution 
equations is a challenging task. The main purpose of this paper is to report a new, simple, 
and efficient numerical procedure for this model. 

Since the corresponding initial value problem is typically stiff, implicit time stepping 
methods should be implemented. For the Simo & Miehe version of the MF, a discretized 
problem can be obtained using the operator split technique in combination with expo- 
nential mapping and formulation in principal axes as described in [63] and [51] (Eule- 
rian approach). Alternatively, the evolution equation formulated on the reference con- 
figuration can be discretized as described in [11] (Lagrangian approach). In both cases, 
a system of nonlinear algebraic equations is obtained, and a local iterative procedure 
is usually implemented to resolve the resulting nonlinear problem (see, among others, 
[51,44,12,13,49,19,30,55,66,15,22,50,36]). Obviously, such iterative procedures can slow 
down the entire FEM simulation. This problem may become especially important if glob- 
ally explicit FEM is consideredlU This publication is dealing with first-order accurate 
methods only. For the discussion concerning the application of higher order methods, the 
reader is referred to [52,14,5,6]. 

In order to speed up the FEM computations, much attention was paid to the construction 
of closed form solutions for implicit schemes. For instance, a simplified fiow rule under 
the assumption of small elastic strains was considered by Simo and Miehe [62] in order 
to get an explicit update formula for the local implicit time stepping procedure. For the 
same reason, another simplification of the fiow rule in case of small elastic strains was 
considered by Reese and Govondjee [51]. This simplified version was implemented later 
in [26]. Unfortunately, the simplifying assumption of small elastic strains is not valid 
for many materials like plastics, rubber, biological tissues etc. Moreover, if the modified 
Maxwell body is used to capture nonlinear kinematic hardening in metals, a general finite 
strain version of the model must be utilized as wello Another approach to closed form 
solution is based on special assumptions concerning the energy storage. In particular, a 
quadratic logarithmic strain energy (so-called Hencky strain energy) can be assumed in 
order to simplify the numerical treatment of the material model [41] . Unfortunately, this 



^ In the case of explicit FEM, the evaluation of the material routine at each point of Gauss 
integration constitutes the major part of the overall computational effort. 

^ In fact, although the elastic strains in metals are typically small (Se — ?• in Fig. Ic), the 
conservative part (eie in Fig. Ic) of the inelastic strain may become finite. 



assumption would yield unrealistic results in case of large elastic strains. Thus, again, the 
applicability area is limited to moderate elastic strains. In this work, a simple explicit up- 
date formula is presented for the original finite strain version with Neo-Hookean potential. 
Interestingly, this explicit solution for the general case is even more compact and simple 
than the solutions presented in [62] and [51] for the special case of small elastic strains 
or the solution in [41] for quadratic logarithmic strain energy. For the new method, the 
computational effort per single time step is even smaller than the effort required within 
the explicit time stepping. 

The inelastic flow is assumed to be incompressible, and the algorithm presented in this 
work preserves this incompressibility constraint. A classical model of finite strain vis- 
coplasticity which contains the Simo & Miehe version of the MF was considered in [56] . 
As it was shown in [56], the exact solution to the initial value problem is exponentially 
stable with respect to small perturbations of the initial data, if the incompressibility con- 
straint is not violated. For such material models, the numerical schemes which exactly 
preserve the incompressibility are advantageous due to the suppressed error accumulation 
[56]. This theoretical result is confirmed by numerical tests presented in the current paper. 

Dealing with the constitutive equations written in Lagrangian form, it can be shown that 
they are invariant under the isochoric change of the reference configuration [59] . The same 
invariance property can be formulated for the numerical solution as well. Obviously, the 
numerical algorithms which exactly retain this invariance property are advantageous. In 
this work, it is proved that the advocated algorithm retains the invariance of the solution. 

We conclude this introduction with a few words regarding notation. Throughout this 
article, bold-faced symbols denote first- and second-rank tensors in M^. A coordinate- free 
tensor formalism is used in this work [25,55]. In this work, 1 stands for the second-rank 
identity tensor. The deviatoric part of a tensor is defined as A° := A — |tr(A)l, where 
tr(A) stands for the trace. The material time derivative is denoted by dot: ^A = A. The 
overline (■) denotes the unimodular part of a tensor such that A = (det A)^^/^A. The 
inverse of transposed tensor is denoted by A^""". The positive definiteness of a tensor A 
is symbolically denoted by A > 0. 



2 System of constitutive equations 

2. 1 Lagrangian formulation 



Let us consider a finite strain model of Maxwell fluid. This model is covered as a spe- 
cial case by the viscoplasticity model presented by Simo and Miehe [62]. The Lagrangian 
formulation of the model follows the presentation of Lion [38] . We start with the multi- 
plicative decomposition of the deformation gradient F into the elastic part Fe and the 



inelastic part Fi 

F = FeFi. 

Along with the well-known right Cauchy-Green tensor C = F"'"F, we introduce the inelas- 
tic right Cauchy-Green tensor as 

Ci = F^Fi. 

The elastic right Cauchy-Green tensor Cc and the elastic Green tensor Fg are defined by 

Ce:=FjFe, fe = i(Ce-l). 

Next, we introduce the inelastic velocity gradient L; and the covariant Oldroyd derivative 
(with respect to the intermediate configuration) 

The inelastic Almansi strain tensor Fi and the inelastic strain rate tensor Dj are defined 
through 

After some straightforward computations (cf. [16]), one gets 

A 

Di =ri . (1) 

Let T be the Cauchy stress tensor. The Kirchhoff stress tensor S, the 2nd Piola-Kirchhoff 
tensor S operating on the intermediate configuration and the 2nd Piola-Kirchhoff tensor 
T operating on the reference configuration are defined through 

S := (detF)T, S := Fe^^SF;^, t := F-^SF-^. (2) 

Let ip be the free energy per unit mass. It is postulated in the Neo-Hookean form as 

PH^(Ce) = ^(tra-3). 

Here, pi^ > stands for the mass density with respect to the reference configuration; 
/x > is the shear modulus. Next, noting that Ce = 1 + 2re, a hyperelastic stress-strain 
relation is introduced on the intermediate configuration: 

^ gV^(l + 2fe) 

S = Pr -^ . (3) 

Isothermal processes are considered in this study. The Clausius-Duhem inequality requires 
that the specific internal dissipation 6i remains non- negative (see [16]) 

(5i := — t : E - ^ > 0, (4) 



where E := 7j(C — 1) stands for the Green strain tensor. Using (3) and taking the isotropy 
of the free energy function into account, this inequahty is reduced to 

A 

pj, = (CeS) :f i> 0. (5) 



An evolution equation is postulated so that inequality (5) holds for arbitrary mechanical 
loadings (cf. [38]) 

ri= ^(CeS)°, (6) 

where i] >0 is a material parameter (Newtonian viscosity) l£j In view of (1), an equivalent 
formulation of this flow rule is given by 

A = ^(CeS)°. (7) 

Note that this flow rule is six dimensional since the rotational part skew (Li) drops out from 

A 

the constitutive relations due to the elastic isotropy. Moreover, since tr(ri) = tr(Di) = 0, 
the inelastic flow described by (6) in incompressible. 

In order to simplify the numerical treatment of the model, the constitutive equations can 
be transformed to the reference configuration. In particular, the free energy takes the form 



/^ 



7/> = ^(ccr') = -^(trccr' 

Using (2) and (3), one gets for the 2nd Piola-Kirchhoff stress tensor 



T = 2pR- 



Ci=const 



dC 
Substituting (8) into this relation, one gets 

t = /iC-^(CCr')°. (9) 

Next, we note that 

tr(CeS) =tr(Ct). 

Thus, the pull-back of the deviatoric part of the Mandel tensor CgS is given by 
F;^(CeS)°Fi = CtCi - ^tr(CeS) C; ^^ (ct)°Ci. 



^ In some works, the viscosity r] is replaced by the parameter x > such that x = l/rj. 
Furthermore, the relaxation time can be introduced by r = r]/^. The case of process-dependent 
viscosity was considered, among others, in [37,33,38,57,31]. 



Using this relation in combination with the evolution equation (6), we get 
Q = 2F;^ f i Fi ^^ lF;^(CeS)°Fi ^^ -(ct)°Ci. 



Taking (9) into account, we have 



Ci = -(CCr^)°Ci. (10) 



The system of constitutive equations (9) and (10) is closed by specifying initial condinions 

Ci\t=to = C?. (11) 

Note that the exact solution to the evolution equation (10) has the following geometric 
property 

Ci(t) em if c° G M, (12) 

where the manifold M is a set of symmetric unimodular tensors 

M := {a G Sym : detA = l}. 

It follows from (12) that Ci remains positive definite if Cf > Ol£j 

Remark 1 Interestingly, an evolution equation in the form (10) was presented for an 
alternative model of the MP by other authors (see equation (10.147) in [16]), although 
the flow rule considered on the intermediate configuration is given by Di = -S, which 
differs essentially from the flow rule (7). Since both evolution equations coincide, the 
approach advocated in this paper can be applied to the version presented in [16] without 
any modifications whatsoever. 

2.2 Eulerian formulation 



Now let us check that the model presented in this section coincides with the model of 
Simo and Miehe [62], formulated within the Eulerian approach. Here, elastic isotropy is 
considered. First, note that 

2Di = FrTCiFri. (13) 

Due to the elastic isotropy, according to the evolution equation (7), Di commutes with 
Ce- In particular, one gets from (7) 



Ce2DiC-i = -(CeS)°. (14) 



^ The eigenvalues of Ci are continuous functions of time. Therefore, if all of the eigenvalues were 
positive at some time instance and their product remains constant in time, then they remain 
positive. 



Substituting (13) into (14) and taking into account that j|(Cj ) = — C; CjC; one gets 

- Fy^^{Cr^)F^F-^F-'F-^ = ^(CeS)°. (15) 

For what follows, we introduce the elastic left Cauchy-Green tensor Be = FeFj, and note 
that 

S° = F-^(CeS)°Fr. 

Multiplying (15) with F^T^ from the left and with Fj from the right one gets using (2.2) 

_Fl(Cr')F-B.-4s°. (16) 



By introducing the covariant Oldroyd rata ^ I 



C/^f A ^ ._ T? /'T?-l A T?-T\TriT _ A T A A T T 



D(A) = i^lA) := F— (F-^AF-^)F^ = A - LA - AL 

at 

one gets 

D(Be) = ^(Be) = F|(Cr^)FT. 

Thus, the flow rule (16) takes the well-known form 

- ^(Be)B;i = -S°. (17) 

This equation was covered as a special case by Simo and Miehe [62] (see equations (2.19a) 
and (2.26) in [62]) as well as by Reese and Govindjee [51]. In case of the Neo-Hookean 

potential (8), we have S = S° = /u(Bc)°, and the evolution equation (17) is reduced to 

- =^(Be)B,-i = ^(Kf. (18) 

V 



3 Time stepping algorithm 

3. 1 Explicit update formula in Lagrangian formulation 



Let us consider a typical time interval (t„,t„+i) with At := t„+i — t„ > 0. By "Ci 
and "^^Ci we denote numerical solutions respectively at t„ and t„+i. Suppose that the 
deformation gradient "+^F at time instance t„+i is known, and "Ci G M is given. The 
unknown ""'"^Ci G M is estimated as the unimodular part of the solution provided by 



In order to sress that this Oldroyd rate is also a Lee derivative, it can be denoted by =-4f- 



the classical Euler-backward method (EBM). In other words, let "■+icf^'^ be the EBM 
solution. After the subsequent correction, the solution is given by 



n+l. 



^1/3 



^Ci:="+iC™^^=(det("+^C™^)) ^""+iC™M_ (^g) 

Let us derive this solution. First, adopting EBM, the discretized version of (10) reads 

where "+^0 = "'+iF^ «+ip jg given. Introducing abbreviation 

P := i^ ^^(n+l^n+l(.EBM)-l)^ (21) 

3 r] 



(20) can be rewritten as 



n+l/-iEBM _ n/-i , f^^n+1^ _ n n+l/iEBM 



n+l/-^EBM _ -'- fnr^ , ^^I^n+IT^ 



Thus, obviously 

(^EBM _ ^^ /n(^. ^ ;^:2:«+lC ) . (22) 

1 + p ^ r] ^ 

Finally, noting that YXgA = A for any A, and substituting (22) into (19), the following 
explicit update formula is obtained 



„+!(.. 




= nQ. + ^^/^n+lC 



(23) 



Remark 2 Note that the solution ""'"^Ci is obtained directly, without computing "+icp^^. 
On the other hand, a simple explicit update formula can be derived for EBM as well (see 
Appendix A). 



3.2 Properties of the algorithm in Lagrangian formulation 



Recall that the exact solution to the evolution equation (10) has under proper initial 
conditions the geometric property Ci G M. The numerical algorithms which exactly pre- 
serve this property are referred to as geometric integrators. It was shown in [56] that 
such integrators allow to suppress the error accumulation, which is especially important 
for the simulation of "long" processes. Therefore, geometric integrators are advantageous. 
Obviously, the geometric property (12) is exactly retained by the numerical solution (23). 

Remark 3 Observe that the method (19) corresponds to EBM with a subsequent cor- 
rection. Another modification of the classical EBM was considered by Helm [19] in order 



to obtain a geometric integrator. Furthermore, in the paper by Vladimirov et al. [66], 
two other modifications were considered and the plastic incompressibihty was enforced 
by introducing the additional equation det Cj = 1 with an additional unknown scalar 
variable. In contrast to the explicit update formula (23), a local iterative procedure was 
implemented in [19] and [66]. 

Another property of the algorithm is as follows. As already mentioned in Section 2.1, the 
exact solution Ci must be positive definite. Since "Ci and ""'"^C are positive definite, the 
numerical solution "+^Ci given by (23) is positive definite as well. Indeed, the sum of 
two positive definite tensors is positive definite, and the projection operator (■) retains 

this property. Furthermore, due to the positive definiteness, det("Ci H ^"+^0) > 0. 

Thus, the unimodular part of ("Cj H — ti£n+iQ^ ^^ ^23) is well defined for all At > and 
"'''^Ci it is a smooth function of At. For At > 0, the solution "+^Ci ranges smoothly from 
"Ci to "+^C. Thus, there is no danger of "wild oscillations" of the solution, which are 
typical for explicit time stepping schemes dealing with large time steps. The algorithm is 
unconditionally stable since the solution remains finite for arbitrary time steps. Like the 
classical EBM, the algorithm is first order accurate (see Appendix B). 

In case of a relaxation process with C = const, the Clausius-Duhem inequality (4) re- 
quires that the free energy is a decreasing function of time. Let us analyze the dissipation 
properties of the presented algorithm in case of stress relaxation. We consider a local 
relaxation process with C = const and Ci = "Cj as prescribed initial condition at t„. In 
that case, it can be shown that ip{C ""*"^Ci ) is a monotonically decreasing function of 
At (see Appendix C). Thus, the relaxation property of the exact solution is qualitatively 
reproduced by the numerical solution. 

Note that the evolution equations (10) are invariant under isochoric change of the reference 
configuration (see [56,59]). More precisely, let Fq be a constant tensor such that det Fq = 1. 
The new reference configuration is introduced as K^'^'" = FqK. Thus, the corresponding 
deformation gradient (relative deformation gradient) is then given by F°^™ = FFq"^. Along 
with the "old" quantities C and Cj, we consider their new counterparts 

In that case, it can be easily shown that the evolution of Cf^"^ is governed by (10) if 
all "old" variables are formally replaced by their new counterparts. The same invariance 
requirement also can be formulated for the time-stepping procedure: For a fixed time step 
At, we denote the numerical solution by "^^Ci = 9T("^^C, "Ci). The numerical scheme 
9T is called invariant under the reference change if 

g^^p-T n+1^ p-1^ p-T n(^. p-1^ _ ¥^^^^+^0, "Ci)Fo"^ (24) 

This property can be checked for the scheme (23). Indeed, since det Fq = 1, we have 



Fo-^"CiFo-^ + ^ Fq-^ "+1C Fq-^ = F,^ -Ci + ^ "+1C F^\ 
1] 1] 
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Remark 4 One possible application of the invariance property (24) is as follows. By 
choosing Fq := ("Ci)"*^/^ we get "C°®™ = 1. Thus, any algorithm computing D^("+^C, 1) 
would be sufficient to restore the general scheme 9T. In other words, it is sufficient to 
perform time-stepping for zero initial conditions for the inelastic strain, which corresponds 
to "Ci = 1. 

Within implicit, deformation driven finite element procedures, the consistent tangent op- 
erator is required for iterative solving for global balance of linear momentum [64]. An 
explicit expression for the consistent tangent is presented in Appendix D, and its symme- 
try is proved. 

It can be shown that the numerical solution (23) corresponds to the exact solution of a 
relaxation problem with a fixed C at a certain time instance t ~ tn+i- More precisely, 
consider the initial value problem (10), (11) with C(t) = const. Then the exact solution 
can be represented as Ci(t) = C" + ip{t)C with ip{t^) = and ip{t) ^ ^^^^^^ as t ^ t° 
(see Appendix E). 



3.3 Explicit update formula in Eulerian formulation 



In this subsection we rewrite the explicit update formula (23) in terms of tensors which 
operate on the current configuration. Toward that end, for each time interval (tri,^n+i), 
we introduce the so-called trial elastic left Cauchy-Green tensor under the assumption of 
a frozen inelastic flow. Having in mind that Be = F C;"^ F""", we get the trial elastic strain 
in the form 

n+l-ptrial n+lrp n/— ( — 1 n+lrpT 

Multiplying (23) with "+iF from the left and with "'+^F from the right, we get 



"+lg _ ('n+lgtriaU-l i ^^ ]_ 



Taking the inelastic incompressibility into account, one gets 



n+lg-l 




\ 1 1 1 


- (det"+^F)-2/3(n+im-i)-i + ""'^ 1 



(25) 



This is an explicit update formula for the evolution equation (18). Note that, similar to 
the product formula consistent with the operator split considered by Simo [63], """"^Bc is 
co-axial with "+13*''''^^ Within a time step, the explicit update formula (25) predicts the 
same stress response as the formula (23). 
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3-4 Exponential mapping: Lagrangian and Eulerian formulations 



Alternatively to the explicit update formula considered in this paper, a well-known implicit 
method based on the exponential mapping can be adopted: 



n+l 



Ci = exp 



^if' /n+lp n+lp-lND 



^Q. 



(26) 



By some algebraic computations, it can be shown that the invariance relation (24) holds 
for this method as well. Next, a proof that the symmetry property is retained by the 
algorithm even in case of a more general material model is presented in [54]. Moreover, 
the incompressibility of the inelastic flow is exactly retained. Thus, for the corresponding 
numerical solution, "'^-'^Ci G M. 

Observe that the numerical scheme obtained using the operator split technique (cf. Simo 
[63]) can be derived directly from (26). Indeed, by inverting both sides of (26), we get 



n+lf-i — l n/~( — 1 



Cr' = "Cr' exp 



^H^/n+l-^ n+l^-l\D 



(27) 



Since the inelastic incompressibility is exactly retained by (26), det(""'"^Bj, ^) = (det "+iF)^. 
Therefore 



n+lrp-T /n+17=^ n+l/-i-l\D n+l-rri-T _ m+lfi \D 



L/~( n+i/~( 



B, 



(2^ 



Multiplying (27) with "+iF from the left and with "+iF"'" from the right, we get using 
(28) 



n+l-D n+1-Dtrial 



^B 



B^"^^ exp 



^("+iBe)° 
V 



Therefore, ""'"^Be commutes with "+13*''"^^ Thus, this equation can be rewritten in the 
form 



n+l 



Be = exp 



^r-^^I. 



n+l-ptrial 



(29) 



The well-known implicit update formula is restored (see, for instance, equation (44) in 
[51]). We stress that the nonlinear equations (26) and (29) represent one and the same 
method, written in two different formulations. For a given time step, these equations are 
equivalent since they predict the same stress response. 
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4 Numerical results 



4.1 Accuracy testing for a single Maxwell element 



In order to test the accuracy of the explicit update formula (23), we consider a local 
loading program in the time interval t G [0, 300] (time is measured in seconds) 



F(t) = F'(t), (30) 

where F'(t) is a piecewise linear function of time t such that F'(0) = Fi, F'(IOO) = F2, 
F'(200) = F3, and F'(300) = F4 with 

Fi := 1, F2 := 2ei (g) ei + ^(ea ® eg + eg (g) eg), 

v2 

F3 := 1 + ei (g) e2, F4 := 2e2 (g e2 + ^=(ei (g ei + eg (g eg). 



V2' 



More precisely, we put here 



- t/100)Fi + (t/100)F2 if t e [0, 100] 
F'(t) := <( (2 - t/100)F2 + (t/100 - 1)F3 if t G (100, 200] . 

- t/100)Fg + (t/100 - 2)F4 if t G (200, 300] 

The reference configuration is assumed to be stress free at t = 0. Thus, we put Ci|(=o = 1. 
The following values of the material parameters are used: f] = 400 MPa s, /i = 40 MPa. 
The numerical solution of the initial value problem obtained with extremely small time 
step (At = 0.001) will be referred to as the exact solution C?^"'^*. The numerical solutions 
with At = 1 and At = 0.5 are denoted by Cp'""''. The error ||Cr™"'' - Cf^"'^*!! is plotted 
versus time in Fig. 2 for three different methods. The explicit update (23) is abbreviated as 
Euler-backward method with subsequent correction (EBMSC). Additionally to EBMSC, 
the classical EBM method and the exponential method (EM) are considered in this sub- 
section. Since these three methods are first order accurate, the error is proportional to 
At. Moreover, in accordance with [56] , there is no error accumulation in case of geometric 
integrators (EBMSC and EM). More precisely, the error is uniformly bounded by CAt 
where the constant C does not depend on the size of the entire time interval [56]. Next, 
since the incompressibility condition is violated by EBM, the geometric property (12) 
is not preserved and the numerical error tends to accumulate over time. The geometric 
integrators EBMSC and EM are equivalent in terms of accuracy. The simulated stress 
response is presented in Fig. 3 for three different time steps At. 

Next, the accuracy is tested for varying viscosity rj with a fixed /i = 40 MPa, At = 10 and 
loading programm given by (30). The largest relative error is observed for the smallest 
values of 77, although the overall stress level tends to zero (see the top part of Fig. 4). 
Indeed, for smaller 77 the problems becomes stiffer. Although the implicit schemes are 
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Fig. 2. Plots of the numerical error pertaining to three different algorithms: classical Euler 
backward (EBM), exponential method (EM), and the Euler backward with subsequent correction 
(EBMSC). Explicit update formula for EBMSC is given by (23). All methods are first order 
accurate. Geometric integrators (EBMSC and EM) prevent the error accumulation. 

robust, smaller time steps are still required to reduce the numerical error. In case of large 
rj, the inelastic flow is insignificant and the material response becomes nearly hyperelastic 
(see the bottom part of Fig. 4). 



4-2 Application to visoelasticity within the FEM 



4-2.1 Material model of finite strain viscoelasticity 

In this section, a practical application of the proposed integration algorithm to the mod- 
eling of viscoelastic material response is discussed. As a typical example, a constitutive 
model proposed by Reese and Govindjee [51,10] is considered. The rheological interpreta- 
tion of the model consists of a single hyperelastic spring (Hooke body) for the represen- 
tation of equilibrium stresses and N Maxwell bodies which are connected in parallel to 
the Hooke body to represent viscous effects, see Fig. lb. Each of these Maxwell bodies is 
governed by equations presented in Section 2. The total free energy is given as a sum of 
isotropic functions as follows (cf. [51,10]): 

N 



i, = ^eq(C) + E ^ov,^(CCr^: 



m=l 



where the equilibrium part is modeled by the hyperelastic Yeoh model [68] in combination 
with an additional volumetric part 



p,7Aeq(C) = J2 c.o(trC - 3)" + ^( (detF)i/3 _ 1 ) . 



n=l 



(31) 



Here, Cio, C20, and C30 are material parameters of the Yeoh material; k stands for the bulk 
modulus. Observe that the volumetric stress response depends on the volume ratio, given 
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Fig. 3. Simulated stress response in terms of Cauchy stresses for At = 10 (top), Ai = 5 (middle), 
and At = 1 (bottom). Since both methods are first order accurate, the error is proportional to 
At. 



by det F|_i| Next, the free energy for the m Maxwell body reads as 



PKWow,m — PR'?/'ov,m(CCi;j) — -— (^trCCj„j — 3 



m = l,2,...,iV. 



Here, /i^ > is the shear modulus of the m*^ element and Ci^m is the corresponding 
inelastic tensor of right Cauchy Green type. The evolution of each of these variables is 



^ The volumetric ansatz presented in (31) is implemented in MSC.MARC [42]. Obviously, any 
alternative ansatz for the volumetric part can be used as well. 
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Fig. 4. Simulated stress response in terms of Cauchy stresses for different values of ??: 
7] = 4 ■ lO^MPa s (top), r] = 4 ■ lO^MPa s (middle), and ?? = 4 • lO^MPa s (bottom). The 
relative error is higher for stiffer problems with small r/. 

governed by equations of type (10). The resulting system of constitutive equations is as 
follows: 

N 

J- J- oq ~r / ^ J-ov,™,; 

m=l 

teq= f2cio+4c2o(trC-3)+6c3o(trC-3)^) C°C-^+3A;(detF)^/3 ((det F)i/=^-l) C-\ 



/ -, X D 



/ -, \D 



th 



ra, m = l,2,...,N, 



where rjm > denotes the viscosity of the m Maxwell body. 



16 




1 1.2 1.4 1.6 1 



fe 




-Is^lO 

CL, 



500 1000 

time t [s] 



2000 



viscoelastic model 
equilibrium part 




1 1.2 1.4 1.6 1.8 



2.2 2.4 
F 

■■- T.I 



Fig. 5. Simulated stress response under uniaxial cyclic loading for different strain rates: for 
|F^^| = 1.5 s-i (left) and for |F^.^.| = 0.015 s'^ (right). 

The material model is implemented into commercial finite element software MSCMARC®, 
making use of the user subroutine HYPELA2. Since the total strain is prescribed at each 
point of Gauss integration, the numerical treatment of each of the Maxwell elements is 
fully independent from other Maxwell elements. For each of them, the corresponding evo- 
lution equations are solved with the help of the explicit update formula (23) and the 
consistent tangent operator is computed explicitly (cf. Appendix D). In order to avoid 
volume locking effects, a mixed u-p- formulation is adopted [65] . For the subsequent simu- 
lations, a system of 4 Maxwell bodies is implemented with the material parameters from 
Tab. 1. 

Table 1 

Material parameters of a viscoelastic material (m = 1,2,3,4) 



cio [MPa] 


C20 [MPa] 


C30 [MPa] 


k [MPa] 


N[-] 


/im [MPa] 


r/„ [MPa s] 


0.45 


-0.048 


0.011 


1000 


4 


0.2 


2 • 10™-3 



First, in order to illustrate the stress response of the model, a uniaxial cyclic loading is 
simulated. The loading axis is fixed and coincides with the x-axis. In Fig. 5 the stress 
response under a strain driven loading is shown for two different loading rates. It can be 
seen that for small loading rates the stress response converges to the equilibrium path. 
Next, the simulated stress response for uniaxial relaxation test is represented in Fig. 6. 
Note that if the instant stress lies under the equilibrium stress, the relaxation process 
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Fig. 6. Simulation results for uniaxial relaxation tests: loading program (top) and stress response 
(bottom). 

results in increasing stresses. 



4-2.2 FEM solution of a representative boundary value problem 

A representative boundary value problem is solved in this subsection in order to test the 
stability properties of the proposed integration algorithm. Toward that end, the material 
model from the previous subsection is adopted to simulate the so-called rotary dynamics 
experiment used for the analysis of rubber materials. The corresponding experimental 
setup was originally proposed in [9] and applied later for the hfe-time prediction of rubber 
materials in [28,29]. Within this experiment, a shear loading with rotating axes is applied 
to two cylindrical specimens as shown in Fig. 7i^ 

The used finite element model is adopted from [21]. It represents a single cylinder of 
height /i = 40 mm and of diameter c? = 10 mm (Fig. 8). Eight-node isoparametric three- 
dimensional brick elements with trilinear interpolation and one extra node with a single 
degree of freedom for pressure (MSCMARC element type 84) are used. The mesh consists 
of 7000 elements. Loads are applied only on the upper and lower surfaces in the following 
manner. Both surfaces are rigid and can rotate independently about corresponding rota- 
tion axes. Each rotation axis goes through the center of the corresponding surface in the 
normal direction. For the lower surface, the rotation axis is fixed and the rotation is free. 
For the upper surface, an initial displacement of the surface is submitted in ^/-direction 
as shown in Fig. 8, such that a simple shear loading is applied to the specimen. In the 



The animated version of the experiment can be seen at http://youtu.be/eNgjGE7upYY . 
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Fig. 7. Experimental setup for simple shear loading with rotating axes [28]: initial state (left), 
monotonic shear loading (middle), shear loading with rotating axes (right). 
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■^ ^ 
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Fig. 8. FEM model of the rotary dynamics testing: initial configuration (left) and deformed 
configuration (right). 

next loading step, the corresponding axis is fixed and a rotation of the upper surface is 
prescribed with a constant angular velocity ipz = 0.27r/s. Thus, the sample is subjected 
to a non-proportional cyclic loading. Finally, after 300 s, the rotation is stopped. 

The numerical simulation was performed with the time step At = 0.25 s. The material 
parameters of the viscoelastic material are taken from Tab. 1. No convergence difficulties 
were observed during the simulation whatsoever. The simulated reaction forces F^ and Fy 
are shown in Fig. 9. After approximately 10 revolutions of the sample, nearly stationary 
solution is observed. Finally, during the relaxation step, the reaction force F^ disappears 
and the force Fy tends to a certain equilibrium value. 



5 Discussion and conclusion 



A closed from solution for fully implicit integration algorithm is proposed. The algorithm 
corresponds to the classical Euler-backward method with a subsequent correction to en- 
force the incompressibility of the inelastic flow. The scheme is highly efficient, since no 
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Fig. 9. Simulation results for reaction forces. FEM simulation performed in MSC.MARC using 
explicit update formula (23). 

local iteration procedure is required. The numerical solution is well defined even for large 
time steps At. The resulting integration algorithm is unconditionally stable and first order 
accurate. The algorithm shows a similar accuracy as the well-known exponential scheme 
(EM). The consistent tangent operator is symmetric. The following properties of the exact 
solution are retained by the algorithm (23): 

i symmetry and incompressibility: "+^Ci G Sym, det("+^Ci) = 1, 
ii positive definiteness: "+^Ci > 0, 
iii the solution "+^Ci is a smooth function of At, 
iv for the stress relaxation with C = const, the free energy ip{C "^^Ci ) is a monotoni- 

cally decreasing function of At, 
V numerical solution remains invariant under the isochoric change of the reference con- 
figuration. 

No error accumulation is observed due to the exact preservation of the incompressibility 
constraint, in accordance to the theoretical results from [56]. 

The explicit update formula is derived for a special variant of the finite strain Maxwell 
fluid which is widely adopted in material modelling. Some modifications of the algorithm 
are possible to cover more general material behavior. In particular, the case of process- 
dependent viscosity can be considered. Moreover, the Neo-Hookean type of hyperelasticity 
can be replaced by more general constitutive assumptions. Application of the method to 
elasto-plasticity with different types of nonlinear hardening seems promising. 
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Appendix A (explicit update formula for EBM) 



Starting from (22), we derive an explicit update formula for the classical Euler-backward 
method. It is sufficient to obtain a closed-form relation for f3 which appears in (22). First, 
we abbreviate: 

V 
By computing the inverse of both sides of (22), we get 

(^-+^C™)-' = {1 + /3)^-\ 

Substituting this result into the definition (21), we obtain a linear equation with respect 
to/3 



/3 = -^^tr (l + /3)"+^C$ 
Si] ^ 



After resolving it, we get 

/3 



1^ tr("+iC#-i) 
1-ii^ tr(«+iC#-i)' 

3 r; V / 



Finally, substituting this into (22), the EBM solution is given by 



n+l/-iEBM _ /i I fjAJ: . fn+lr^^-l 



3 r] 



Appendix B (convergence rate) 



Let us show that the method (19) is first order accurate. Toward that end we consider a 
time interval (t„,t„_|_i). Let Ci(t„+i) be the exact solution to the problem (10) with the 
initial condition Ci\t=t = "Ci. It is sufficient to show that there exists C < oo such that 



in+l 



Ci - Ci(t„+i) II < C {Atf as At -^ 0. 



It is well known that the EBM is first order accurate. Moreover, there exists Ci < oo 
such that 

||„+i(.EBM _ Ci(t„+i)|| < Ci (At)2 as At ^ 0. 

On the other hand, (det(X))^^/'^ is a smooth function of X in vicinity of the exact solution 
Ci(t„+i). Therefore, for sufficiently small time steps, there exists a constant C3 < oo such 
that 

|(det("+iC™^)-i/3 _ i| = |(det("+iC™^)"i/3 _ det(Ci(t„+i))-i/=^| < 

Cs f+^C™ - Ci(t„+i)|| < C1C3 (At)2 as At ^ 0. 
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Thus, for small At, there exists C4 < oo such that 



||n+lcEBM _n+l(.EBM|| ^ \(^aetC+'C™)-'/^ - 1| \\r^+^C™\\ < C4 (Atf. 

Finally, using the triangle inequality, we get the required estimation 



in+l 



Ci - Ci(t„+i)ll = r+'C™ - Q(Wi)|| < 

||n+lcEBM_n+l(.EBM|| ^ ||n+l(.EBM _ (..(^^^^) || < (^7^ + Ci) (At)^ 



Appendix C (energy release during relaxation) 

Let us show that ip{C "^^Ci ) is a monotonically decreasing function of At for a fixed 
C, where ""'"^Ci is given by the explicit update formula (23). We make use of the fact 
that the free energy ^{C Ci~^) is invariant under the isochoric change of the reference 
configuration. Thus, due to the invariance of the integration algorithm with respect to 
the reference change, we may assume that C = 1. In that case, the update formula takes 
the simple form 

"+^Ci = #, $:="Ci + ^l. (32) 

V 

For C = 1, it follows from (8) that 



c^At^^ ' 2pa dAt^ ^ ' 'J 2pR V dAV ' 



(33) 



Note that 



d fi d 



-$ = ^1 



rfAt 1] ' rfAt^ ' 3r] 

Thus, differentiating (32) i and using (34), we get 



(det #)-^/=' = -i^(det #)-^/^ tr(#-i). (34) 



d n+l 



Q = (det$)-^/=^^[l - -tr($-i)$ 



dAt V 

Substituting this result into (33) and taking into account that "-^^C^^ = (det #)^/^#~^, 
we get 

7 2 1 

— -^(C "+iCi~') = -(det$)i/3 -^ tr($-i[l - -tr(#-^)#]#-i) = 



c^At"^' ' ' ' 2p^7] ' ' 3 

_(det#)i/3 _^tr(#-i(#-i)D) = -(det$)i/3 -^||(#-i)°f < 0. 

In particular, the free energy of the numerical solution for relaxation processes is a de- 
creasing function of At. 
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Appendix D (consistent tangent operator) 

Combining (9) and (23), we define a function Ti(C) tlirougli 



^Ur^U^ . Mt^\-^ 



D 



Ti(C)=^C-^|C("Ci + ^^C) I . (35) 



Next, we abbreviate 



$:=-Q + ^C, Q(C):="Ci + ^C = #. (36) 

rj 7] 

Using notations from [55], the consistent tangent operator is explicitly given by 



at.,c),.,a(.c-Mccrr,_^^_.^^_.(^^_.^, 



dC dC 

,.(detC)-"=C-'.|--(CCri)i' »c-' + P : (l ■ Cr' + C • i^)|, (37) 

^.-(cr'0cr'):§. (38) 

§'SNdet*)-'/»(l-i(# «*-)):» (39) 

»«S'^(detC)-A.(l-l(C«C-')). (40) 

Now let us prove the symmetry of the tangent operator — ^^ ' . First, note that the ma- 
terial tangent resulting for frozen inelastic flow (Ci = const) corresponds to the symmetric 
tangent operator of a hyperelastic material: 

dUc-^ (CCri)D) 
oC 

(9T (C) 
Thus, the symmetry of — ilX ' is equivalent to the symmetry of T defined through 



9ti(c) 5(/iC-i (ccr^ 

\Ci=const 

dcr' 



"- ■ o/^ o/-^ \Ci=const 

oC oC 



/i(detC)-^/-^ C-^-h: (C 



dC 
It follows from (38)-(40) that for any Y G Sym 

^ : Y = -^{^-\YC'rcfcr\ (41) 
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Note that tr(AB) = tr(BA) and tr(AB°) = tr(A°B) = tr(A°B°) for arbitrary A, B. 
Using these properties, we get for any X, Y G Sym 

/.(detC)-V3tr|(XC-^)°(c ^ : y) I ^^^ 
_ /^ (detC)-i/3 tr|(XC-i)°('c(#-i(YC-i)°C)°CriU 
_ ^^!^ (detC)-i/3 tr|cr^(XC-i)Dc(#-i (YC-1)DC^° 



_ ^f!^ (detC)-2/3 (det#)i/3 trj (#-i(XC-i)°C)° (^-^ (YC-i)°C)°|. 

Obviously, this expression is symmetric with respect to X and Y. Thus, X : T : Y = Y 
T : X. The symmetry of the consistent tangent — }^ ' is thus proved. 



Appendix E (trajectory of the exact solution) 



Consider the initial value problem (10), (11) with C(t) = const. Suppose that the numer- 
ical solution at time instance t„ is given by 



"Ci = C? + "^C (42) 

with some suitable "(/? > . Substituting this into the update formula (23) we get 



/^^^ ^ ^n , ^ ^ , ^^^^ 



«+iCi = «Ci + ^^^ C = C? + "v'C + ^^^C = C? + "+VC 

with a suitable "+-'^(/9 > 0. Since the assumption (42) is satisfied for n = 1, it is satisfied 
for any n = 1,2,3, ... . Finally, it is known that the numerical solution converges to the 
exact solution as At — )■ 0. Therefore, for any t > t°, we get from (42) 



Q(t) = C? + v^(t)C. (43) 

Thus, the trajectory of the exact solution is given by (43). By substituting this relation 
into (10), the following initial value problem for ip is obtained 

vi = ^(det(C° + v^C))'^', v\t=to = 0. 
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In other words, the parametrization (43) allows to reduce the six-dimensional flow rule to 
a one-dimensional one. 

Finally, we have (/? — )> as t — )■ t°. Thus, since det(C[') = 1 we get det(C° + ip{t)C) ~ 1 
as t — 7> t°. Substituting this into (5), we get ip(t) ^ ^^ " ' as t — ?> t''. 
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